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We present a detailed study of the non-Markovian two-state system dynamics for the regime of 
incoherent quantum tunneling. Using perturbation theory in the system tunneling amplitude A, and 
in the limit of strong system-bath coupling, we determine the short time evolution of the reduced 
density matrix and thereby find a general equation of motion for the non-Markovian evolution 
at longer times. We relate the nonlocality in time due to the non-Markovian effects with the 
environment characteristic response time. In addition, we study the incoherent evolution of a 
system with a double-well potential, where each well consists several quantized energy levels. We 
determine the crossover temperature to a regime where many energy levels in the wells participate 
in the tunneling process, and observe that the required temperature can be much smaller than the 
one associated with the system plasma frequency. We also discuss experimental implications of our 
theoretical analysis. 



I. INTRODUCTION 

It is difficult to overemphasize the importance of the 
dissipative dynamics of a two-state system (TSS). In gen- 
eral, standing as a first hand approximation of a much 
rather complex level structure, the model of a TSS cou- 
pled to a dissipative environment 1,2 has been successfully 
applied to several physical systems. Indeed, the dissipa- 
tive TSS dynamics is the paradigm for the study of super- 
conducting devices containing Josephson junctions j 3 - two- 
level atoms in optical cavities^ electron transfer in bio- 
logical and chemical systems^ and semiconductor quan- 
tum dots^ to name just a few. 

Despite its simplicity, the description of the TSS dis- 
sipative dynamics imposes great theoretical challenges, 
especially when considering non-Markovian processes. 
This is the case for the analysis of the environment low- 
frequency noise spectrum, since the long-lived feature 
of its fluctuations does not allow for a "memoryless" 
bath (Markov) approximation. In the context of a weak 
TSS-bath coupling, theoretical efforts have been made 
to quantify the low-frequency effect for both spin-boson^ 
and 1/f noised models. 

Furthermore, it has been largely demonstrated that 
low-frequency noise plays important role in the deco- 
herence process of superconducting devices containing 
Jospehson junctions. 10 ^ 1 ' 12 ! 13 ' 14 Since those devices are 
seen as promising candidates to the physical implemen- 
tation of a quantum bit, this subject has rapidly grown 
in interest and several studies on describing the micro- 
scopic origin and characterizing the low-frequency noise 
in such devices have already been put forwardj 15 > 16 ' 17 i 18 

Understanding the evolution of a TSS also plays an 
important role in understanding the performance of 
an adiabatic quantum computer— in the presence of 
noisei 20 ' 21 ! 22 ' 23 This is because for many hard problems 
the bottleneck of the computation is passing through a 
point where the gap between the ground state and first 
excited state is very small. Near such an energy anticross- 
ing, the Hamiltonian of the system can be truncated into 
a two-state Hamiltonian^ 3 - and in the regime of strong 



coupling to the environment the two-state results dis- 
cussed in this paper may be directly applied. 

Here, following a previous work^ we put forward a 
detailed study of the TSS dissipative dynamics in the 
presence of low-frequency noise, for the regime of strong 
TSS-bath coupling. We show that, for the regime of small 
tunneling amplitude A, dephasing takes place much ear- 
lier in the evolution, leading the system to incoherent 
quantum dynamics. We employ such a property to de- 
rive equations that describe the non-Markovian evolution 
of the system's density matrix. 

The paper is organized as follows. In section [H] we 
present the system Hamiltonian and a formal solution 
for the time evolution operator. Assuming second or- 
der perturbation theory in A, we calculate in section 11111 
the short-time dynamics of the system reduced density 
matrix elements. Section IIV1 presents a discussion re- 
garding the non-Markovian behavior of the system when 
the environment is in equilibrium. We determine condi- 
tions under which the system reaches the detailed balance 
regime. Section [V] provides a systematic derivation of an 
equation of motion for the system evolution, which in 
general is non-local in time. We also discuss regimes in 
which the equations governing the diagonal part of the 
density matrix become t-local. Considering a double-well 
potential, section I VII puts forward the analysis of intra- 
and interwell transitions and situations where a classical 
mixture of states participate in the quantum tunneling 
process. Finally, section Ivni presents our concluding re- 
marks. 



II. SYSTEM HAMILTONIAN 

We start by considering an open two-state system with 
Hamiltonian 

H = ~[A(t)a x + e(t)cr z }-±<j z Q + H B , (1) 

where Q is an operator acting on the environment de- 
scribed by the Hamiltonian Hb- 
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In order to determine the system evolution operator 
U(t2, ti), we proceed through two simple steps. First, we 
write the state vector of the system Hamiltonian (JT]) as 



\m) = 



>H B t 



\ip(t)). (h = ks = 1, through this paper.) 



Thus, one finds that the state vector \<p(t)) evolves in 
41 



time according to i-S-\ip(t)) = [H (t) + V(i)]\ip(t)), where 



ff (i)=-^K-^Q(t), V(t) = ~A(t)a x , (2) 



and Q(t) 



' J ' Bt Qe lHst . The environment is assumed 



to feature fluctuations following Gaussian distribution, 
therefore all averages can be expressed in terms of the 
correlation function or its Fourier transform, the spectral 
density: 



S{u>) = / dt e^(Q(t)Q(0)), 



(3) 



hence we do not need to specify Hb*^ 

The next step is to make use of the interaction pic- 
ture, considering V(t) as the perturbation. The state 
vector in the interaction picture is defined by \ipi(t)} = 
C/q (t)\(p(t)}, and any operator O is transformed by 

di{t) = ul{t)du Q {t), with 
u (t) = Te -*fi H °vw 



= Texp^ J\e(t') + Q(f)}dt' 



(4) 



where T denotes the time ordering operator. Now, the 
state evolution is determined by the interaction potential 



ffj(t) = ~A(t)ff x (t), 



(5) 



where a x (t) — UQ(t)a x Uo(t). The time evolution opera- 
tor in the interaction representation reads 



_ -i f* 2 HAt)dt 

U I (t 2 ,t 1 )=Te J 'i 



(6) 



Finally we can write a formal solution for the complete 
time evolution operator as 



tf(t 2 ,ti) = Te 



-i f' 2 H(t)dt 



e -'"^U^U^MWlihy 11 ^. (7) 



In this paper, we are interested in the strong coupling 
regime in which the r.m.s. value of the noise 



2^ 



S(u) 



1/2 



(8) 



is much larger than the tunneling amplitude: W 3> A. 
Physically, W is basically the uncertainty in the energy 
bias e(t) and therefore represents the broadening of the 
energy levels. In the above regime, consequently, the 



broadening of the energy levels is larger than the mini- 
mum gap and therefore the gap will not be well-defined. 
On the other hand, as we shall see, for the case of low 
frequency noise, W represents the dephasing rate of the 
system. Thus, the above regime is a limit in which the 
qubit loses quantum coherence before it can tunnel, i.e., 
the dynamics is incoherent. 



III. DENSITY MATRIX CALCULATION 

We would like to study the evolution of the reduced 
density matrix. Let psB{t) denote the total density ma- 
trix of the system plus bath. We therefore have 

PSB(t) = U(t,0)p SB (p)rf{t,0) (9) 
= e~ iHB 1 U a (t) Uj (t , 0)p SB (0) Uj (t , 0) Ul (t)e lH * 4 . 

The system reduced density matrix is defined by p(t) = 
Ti"_b[psb(^)]j where Ttb[...] means averaging over all en- 
vironmental modes. We assume that the density matrix 
at t = is separable, i.e., psb(0) = p(0) <8 Pb, where 



Pb = e 



_ „-H B /T 



is the density matrix of the environment, 
which we assume to be in equilibrium at temperature 
T. Under the assumption of separability of the initial 
density matrix, we consider that the system evolution 
immediately follows an initialization in a definite state, 
implemented, e.g., through a state measurement. 

If A is the smallest energy scale in the problem, we 
can approximate Ui(t,0) by performing a perturbation 
expansion in A, which up to second order reads 



iMt,o)«i + -y o dt'A(t')a x (t') 



4 Jo Jo 



dt'dt"A(t')A(t")a x (t')a x (t"). (10) 



If the time interval t is not small enough to make the 
above integrals small, i.e., t JSJ 1/A, the higher order 
terms in A must be considered in the expansion. 



A. Off-diagonal elements of p 



To zeroth order in A, we have Uj(t,0) = 1, therefore 



PSfl(t) = ^ lHBt U o {t)p SB {0)Ul{t)e 



fill 



For this case, since [Uo(t),a z ] — 0, the a z populations 
of the system reduced density matrix p are constants of 
motion. Therefore, in the representation of the eigen- 
states of a z , a z |0> = - |0)(cr z |1) = |1)), only the off- 
diagonal elements of p present dynamics, which, due to 
the coupling to environment, decay in time. This case 
constitutes the one of a pure dephasing dynamics. It has 
been subject of interest for many areas, where several 
approaches have been used to calculate the off-diagonal 
elements of p. Few examples are the (a) spin-boson 
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model assuming a power-law behavior for the spectral 
density of the bat h 27 i 28 ; (b) spin-fermion mode l 29 ' 30 , and 
(c) spin-two-state fluctuators system 8 -. Here, we consider 
a bosonic bath, but do not have to specify the form of 
the bath spectral density. To quantify this decay, let us 
write the reduced density matrix as 



p(t)= J2 ptf(*)i*>oi- 

i,J=0,l 

We find for the off-diagonal element 



(12) 



Poi(t) = Tr B [(0\U (t)ps B {0M(t)\l)} = poi(O) 

-if\(f)df ^ e - 4 / o t Q(t')dt' Te - 4 / o t Q(t')dt'^ ^ (13) 



x e 



where (...) = Tr b[pb— ] and *T represents the reverse 
time ordering operator. We expand the exponentials, 
group those in the same order, take the average of each 
term, and bring them back to the exponent. Because of 
the Gaussian nature of the environment, it is sufficient 
to expand up to second order in Q. Assuming the envi- 
ronment is in equilibrium, one finds 



j Te~ i /o Q W dt 'Te~%Io Q{t ' )dt ' \ = 

= e _ 3 Sl dt ' f*dt"{Q(t')Q(t")) 



(14) 



Thus, using (jT4|) in (fTBjl . we obtain 
p m (t) =e Jo 



x exp 



dco sm 2 {iot/2) \ 



(15) 



This equation represents a complicated decay rate, which 
is in general not a simple exponential function of t. How- 
ever, in two limits it can be simplified. First, for the case 
of white noise, i.e., S(ui) = 5(0), it gives 



Po 



i(t) = e Jo 2 poi(0). 



(16) 



Which leads to dephasing rate I/T2 = ^S(0), as expected 
for white noise. 

Another interesting limit is when S(lo) is dominated 
by low frequencies so that one can use sin x « x to get 



x[t) = e Jo 2 poi(0), (17) 



where W is the energy level broadening given by ©. 
The decay is now a Gaussian, whose width determines 
the dephasing rate, l/T^ = W. For the case of 1/f noise, 
where the cutoff of S(u>) is not sharp enough, one gets a 
logarithmic correction to the above equation 8 . 



B. Diagonal elements of p 

The evolution of the diagonal part of the density ma- 
trix happens in a time scale much larger than 1/W. The 
complete evolution is given by 

p(t) = Tr B [U (t)U I (t,0)p(0)p B Uj(t,0)U^(t)}. (18) 

Let us assume the initial condition p(0) = |0)(0| and 
try to calculate pu(t). To zeroth order in A, we have 
Pn(t) = as expected, thus we find that the first nonzero 
contribution to pn(t) comes from the first-order term in 

A of HI 



Pll (t) « \ [ dt x ( dt a A(*i)A(t a ) 
4 Jo Jo 

x Tr B [(l|a x (t 1 )|0)p B (0|^(i 2 )|l)] 

= 7 / dti [ dtaA(ti)A(*a) 
4 Jo Jo 

x ^[(lp^u^h^psioiu^Mihm 

= \ fdh f dt 2 k{t 1 )&{t 2 )e^ e(t ' )dt ' 
4 Jo Jo 

tf e -i f* 1 Q(t')dt' Te -± f* 1 Q(t')dt' 



(19) 



where in the second equality we have used a x (t) = 

Ul{t)(T x Ui}{t) = UQ(t)UQ*(t)a x - One can calculate the ex- 
pectation value by expanding the exponentials and keep- 
ing only the the second order terms. The last two lines 
of (TT9"]) become 

1 + \ r ^ r dt " WW)) 

1 Jo Jt 2 
+\ J** dt' £ dt" (Q(t')Q(t")) . (20) 

Substituting the inverse Fourier transformation 
(Q(t')Q(t")) = J ^e- ia, (*'-*")5(w), we find 

;| / ^^)f e M*i-*>) _i +4(8^0^3 -sinorti)] 



2vr to 2 
dio S(u>) 



(coscjt — 1) 



LOT 



2ir to 2 
dto SUo) . . 

jj- f sin lot — 2 sin — cos lot ), (21) 

2n lo z 2 

where r = t 2 — t% and r' = (t% + t 2 )/2. 

If the noise spectral density S(lo) is dominated by low 
frequency noise such that for all relevant modes lot <C 1, 
one can expand the sinwr and coswr in (|2"Tj) to get 

Pll (t) » \ J\t> f dT&{T' + T -W - T -) 
-W 2 r 2 /2^ p (r')r~J:i 2 /2 e(r'+t')dt') 



x e 
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where t = min[2r', 2(t — r')], W is given by {§]), and 

e p (t) = / ^(1-coswt). (23) 

Equation (|22[) conveys the non-locality in time, expected 
for a non-Markovian environment. If within time r ~ 
1/W, e(i) and A(i) do not change much (or even if A(t) 
is a fast but linear exponential function), we can write 
(gH as 

4 Jo J-t 

(24) 

Therefore, for t < 1/A(i), we find the leading term for 
system population rate change given by 



^) /' dr e >[*)-,W]-^- 2 /2. (25) 



If t > 1/W, due to gaussian envelope of the integrand, 
we can extend the integration limits of (|25[) to ±oo, ob- 
taining 

Pii(i)«r p e -[ e (*)-^Wl 2 /2w= 5 (26) 
with the peak value of the functions given by 



^A 2 

w 



(27) 



It is worth recalling that for times t > 1/A, Eq. (jTUJ) does 
not represent a fair approximation to Ui(t,0), hence the 
corrections to Eqs. (|24II26| . due to higher powers of A, 
become appreciable and must be considered. 

In section[V] we present a detailed study for the general 
equation of motion of the reduce density matrix consis- 
tent with (|26|) . However, before we reach that stage, it is 
worth discussing a simpler system with a time indepen- 
dent Hamiltonian, and deriving some general features for 
e p {t) behavior. 



An experimental realization of such a tunneling process 
in a macroscopic quantum device such as a supercon- 
ducting flux qubit is called macroscopic resonant tunnel- 
ing (MKT). The tunneling rates T± are therefore simple 
shifted Gaussian functions described by ([29]) . An imme- 
diate consequence of (1231) is that the shift e p vanishes for 
a classical noise, for which S(u>) is symmetric. Therefore, 
a nonzero value of e p is a signature for quantum nature 
of the noise source. 

If the environmental source is in equilibrium at tem- 
perature T, then the symmetric and antisymmetric (in 
frequency) parts of the noise intensity are related by the 
fluctuation-dissipation theorem: 



S a H = S (w)coth(^) 



(30) 



Therefore the fluctuation-dissipation theorem relates W 
and e p (t) , which are functions of S s and S a respectively. 
Let us first define 



£p0 



duj S(u>) 
2it uj 



(31) 



with V representing principal value integral. In the case 
of low-frequency noise, when all the relevant frequencies 
are small on the scale of temperature T, i.e., ui T, one 
can write coth(u/2T) ~ 2T/u. In that case ©, ([3T?]) . 
and (J3TJ) yield 



W 1 = 2Te 



pO • 



One therefore finds 



e p (t) = e p0 - V J 



dui S(u>) 
2ir oj 



cos Ult. 



(32) 



(33) 



The effect of the last term depends on how small or large 
t is with respect to the time response, ~ UJ c 1 -< °f 
the environment. Here, lo c represents the characteristic 
energy of the environment. To understand this let us 
consider different regimes. 



IV. MACROSCOPIC RESONANT TUNNELING 

Should e p be constant in time and the Hamiltonian ([1]) 
be time independent, one could directly read (|26p as an 
approximation for the equation of motion 



pu(t) w r_poo(i) - T + pn(t), 



(28) 



since the off-diagonal elements of p(t) become negligible 
for times t > 1/W. The rate T_ (r + ) then represents the 
|0) — » |1) (|1) — > |0)) system transition rate. Thus, for the 
regime 1/W < t < 1/A(t), the evolution is described by 
(|2"6")1 . The same argument holds when e p (t) is a function 
of time, but in that case the tunneling rates will be time 
dependent: 



r±(t) 



r p e 



[e±e p (t)] 2 /2W 2 



(29) 



A. Large uj c (short rjj) limit 

If uj c is large compared to 1/t, where t is the typi- 
cal time scale of interest, then the integral in ([33)) cov- 
ers many oscillations of the cosine function and therefore 
vanishes. In that case 



e p ~ e p o 



W_ 

2T ' 



(34) 



consequently being independent of t. For a time indepen- 
dent Hamiltonian, Eq. (|29[) then yields 



r±(e)=r p e-[ £± ^] 2 / 2M/2 , (35) 
in agreement with Ref. [24]- It is easy to see that 



r-(e) 

T + (e) 



(36) 
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which (in the limit A— > 0) is the detailed balance (Ein- 
stein) relation. Therefore, the transition rates ([55]) 
support thermal equilibrium distribution of the system 
states, which is a natural consequence of the fast envi- 
ronmental response. 



B. Small lo c (long tr) limit 

If the environment's response is slow compared to time 
scale of the problem, i.e., lo c -C l/t, the cosine function in 
(|33|) will be close to 1 at all times, making e p sw 0, again 
independent of t. For a time independent Hamiltonian, 
therefore, we get 



r_ =r. 



r„ 



2 /2W 2 



(37) 



Such transitions obviously do not satisfy the detailed bal- 
ance relation and do not lead to equilibrium distribution. 

Indeed, an environment in lu c — > regime behaves as a 
static (classical) noise source. To see this, let us consider 
Hamiltonian |T]) with a static noise source Q that does 
not vary much during the evolution and has a Gaussian 
distribution: 



P(Q) 



2 /2W 2 



2nW 



(38) 



In small A regime, the tunneling rate from state \i) to 
state \j) can be calculated using the Fermi Golden rule 



27r\( i \V\j)\ 2 S(E l ~E J ), 



(39) 



where V = Aa x /2 is the perturbation potential. There- 
fore, for every realization of Q, one finds 



T^(Q) = T + (Q) = ^S(e + Q) 



Averaging over all possibilities of Q, we find 



(40) 



r_ = r i = 



dQP(Q)6(e + Q) 

/2W 2 



(41) 



which is the same as (f3"Tl) . 



C. General lu c (tr) regime 

In general, away from the above two limits, e p (t) is 
a time dependent function given by (|33[) . The explicit 
functionality depends on the spectral density S(u>), espe- 
cially on its characteristic frequency oj c . To see this, let 
us assume a simple spectral density 



5(w) = 



2rjuj 



[1 + (uj/uj c f 



1 



for which analytical solutions is possible. Substituting 
in (j3"3"| . we find 



e p (t) = 



duJ f](l — COSLOt) 



rjuj c 



2tt [1 + (c^K) 2 ] 2 
We can therefore write 



f '(l+tuj c )} 
(43) 



(1 



l (l + tuj c )) 



0. 



UJ c t < 1 
UJ c t > 1 



(44) 



which yields the above two limiting results in the ap- 
propriate regimes with an exponential crossover between 
the two limits. Indeed, the above behavior of e p (t), i.e., 
the crossover from to e p o within time scale ~ l/w c; is 
generic regardless of the functional detail of e p (t). The 
time scale tr ~ 1 /u> c represents the response time of the 
environment to an external perturbation. If t 3> tr, then 
the environment has enough time to enforce equilibrium 
to the system, resulting in 6p — 6po i which is required 
for detailed balance (i.e., equilibrium) condition. On the 
other hand, if t <C tr, the environment cannot respond 
quickly to the system and the equilibrium relation is not 
expected. In that case, we find e p = 0, i.e., the environ- 
ment behaves as a classical noise. In the next section we 
shall see how such behavior results in time-nonlocality of 
the equation of motion. 



V. NON-MARKOVIAN EQUATION OF 
MOTION 

Equation ([26]) gives the short time (1/W < t < 1/A) 
evolution of the diagonal part of the density matrix. As 
soon as t becomes comparable to A, higher order correc- 
tions become important and the second order perturba- 
tion used in Eq. (|26|) becomes insufficient. Instead of 
introducing higher order corrections which is a cumber- 
some task, in this section we take a different path: We 
write a general equation of motion expected for the evolu- 
tion of the density matrix for a system like ours and find 
its parameters in such a way that it agrees with Eq. (|26p 
for short times. 

In general, the equation of motion for the evolution 
of the density matrix is nonlocal in time, reflecting the 
non-Mar kovian nature of the environment. Since the off- 
diagonal elements vanish very quickly (within t ~ 1/W), 
for time scales larger than 1/W, one can write the dy- 
namical equations only in terms of the diagonal elements 
of p. Generally, for a non-Markovian dynamics the equa- 
tion of motion for p(t) depends on the history 



Pii(*)= 



dt'[K-(t, t')p m (t')-K+(t, t') Pll (t% (45) 



where K±(t,t') are nonlocal integration kernels. Let us 
from now onwards consider a time-invariant Hamiltonian 
for which 



(42) pn(t)= 



*'[iMi-<>oo(<')-^+(*-Opii(<')]-(46) 
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If the system starts the evolution from state |0) at time 
t = to, the short time evolution is described by 

pt rt—to 

/>u(t)w / dt'K-(t-t')= / <1tK-(t). (47) 
This should agree with |26|) . therefore 



dTK±(T)=A ± (t-t Q )0{t-t o ). (48) 
where we have defined functions 

A±(*) = r p e -i^ p (t)f/2w\ (49) 

The presence of the ^-function is necessary to ensure 
causality to the system dynamics, since we assume that 
the evolution follows a state initialization at t . Taking 
the derivative of both sides of (|4"8")l . we find 



K±(r) = ^^9(t) + A±(t)$(t). (50) 

Notice that for constant transition rates A±(r) = T±, 
Eq. (JM1 leads to 



pn(t) = r_/9 o(<) - r + pu(t), 



(51) 



which, as expected, is i-local. 

In the limit uj c — > 0, where the change in A± happens 
on a time scale (1/lj c ) much larger than the time evo- 
lution considered here, the time derivative in (|50p can 
be neglected and one obtains ([51]) with transition rates 
r± = A±(Q) = T p e- e2/2w2 , with e p (t) = 0, as expected 
for a static noise. 

On the other hand, in the lo c — > oo limit, variations of 
A± (t) happen in a very short time, hence dA± (r) /dr — > 
for t > tr ~ l/w c . Therefore the i'-integration in (|46|) 
is basically between t— tr and t. If within this short 
range pit') does not change much, one can bring it outside 
the integral. In that case, (|46|) leads to (|5Tj) with r± = 
A±(f -» oo) = T p e-^ ±e ^ 2 / 2W \ with e p (t) = e p0 , which 
is expected in the detailed balance regime. 

Both of the above regimes led to i-local equations for 
the diagonal part of the density matrix. However, for fi- 
nite Lo c , in general, one gets a nonlocal equation in time. 
If the system evolution is slow compared to the time 
scale tr ~ l/w c , one can substitute the Taylor expan- 
sion Pij(t') = pij(t) + it' — t)pij(t) into (|46l) obtaining 



Pll {t) = A_(t)poo(t) - A+(t)pn(t) 



(52) 



where A(i) = A_ (t) + A + (i) . Solving for pn(t), one finds 
(1511) with transition rates 



r± = 



A±(oo) 



A±(oo) 



1- dTTdAir)/di 



•/ °°dr[A(oo)-A(r)]' 
(53) 



where in the last step we have used integration by parts. 
The integral limit was taken to infinity, since the inte- 
grand very quickly vanishes for r > l/u> c . All the nonlo- 
cal behavior is captured in the denominator of (|53p . The 
integrand (|53|) is maximum at r = 0, but very quickly 
vanishes within r ~ l/u; c , hence J °° g?t[A(oo) — A(r)] ~ 
[A(oo) - A(0)]/w c , leading to 



r± 



A±(oo) 



Using 



l-[A(oo)-A(0)]/w c " 
we obtain 



(54) 



A(oo) - A(0) = T P ie~^^ 2 ' 2w2 

= 2T p e-< 2 / 2w2 (e-'W™ 2 cosh-?- - l) . (55) 



Therefore, to the lowest order in T p /lu Ci we get 

2T 



r±(e) « r pe -( £±e ^ 2 /2^| 1 + 



P e -e 2 /2W 2 



( 



e -4o/^ 2 CQsh ^_ 1 
2T 



(56) 



The magnitude and the position of the peak of T_ (e) are 
given by (to the lowest order in T p /uj c ) 



^pcak ~ r_(e p o) ~ Tp(l + r p /u>c), 

2T 



1 + fi£ e -<°' 2W ' 



(57) 



The peak value is enhanced by the nonlocal effects. The 
peak position is also shifted, but by a very small amount 
due to the exponential suppression. Notice that the peak 
becomes asymmetric around its center due to the nonlo- 
cality. 

The nonlocal corrections to the transition rates become 
negligible when T p <C u) c . Also, observe that T p is ap- 
proximately the peak value of the transition rate (|27|) . 
Therefore, nonlocality becomes important only when the 
maximum transition rate F p is of the order of or larger 
than lo c , or equivalently, when the response time (tr) of 
the environment is comparable or longer than the system 
transition time (~ l/T p ). 



VI. MRT IN A DOUBLE- WELL POTENTIAL 

So far we have studied incoherent tunneling in an ide- 
alized two state system. However, for most realistic sys- 
tems, the two state model is only an approximation of 
a more complicated multi-level problem. An example of 
such cases is a system in which the classical potential 
energy has a double-well structure and the kinetic part 
of the Hamiltonian provides quantum tunneling between 
the two wells. Experimental implementation of such 
a system is possible using superconducting flux qubits, 
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which have been studied considerably both theoretically 
and experimentally 3 . 31 i 32 i 33 i 34 . 35 i 36 i 37 i 38 i 39 i 40 . Especially, 
MKT measurements have been performed both between 
ground states as well as excited states of the wellsi 37 ' 38 

In such a double-well system, the energy within each 
well is quantized, with energy level distributions depen- 
dent on the bias energy between the wells. In general, in 
the presence of the environment, a system initialized in 
one of those levels can experience two types of evolution: 
intra- and interwell dynamics. 

The intrawell dynamics are transitions within a single 
well, e.g., when a system excited within a well relaxes to a 
lower energy level in the same well by exchanging energy 
with the environment. Thus, in this case, the system 
dynamics is confined in just one well of the potential, 
with no tunneling to the opposite well. 

It is also possible for the system, depending on the 
tunneling amplitude between the two states, to tunnel 
to an energy level in the opposite well, leading thus to 
an interwell dynamics. If the evolution of the system is 
confined to the ground states of the two wells and it lies 
in the incoherent tunneling regime, then the formalism 
developed herein can describe such an evolution in full 
detail. This, however, is not the only type of evolution 
possible for a double-well system. Here, we also consider 
possibilities that the evolution involves the excited states. 



A. Tunneling between ground states 

At low enough temperatures, the system can only oc- 
cupy the lowest energy states within the wells. In such 
a case, tunneling can occur between those energy levels 
if the levels become in resonance. The probability of the 
system being found in state 1 1) at time t is given by p\\ (t) . 
For a time independent system initialized in state |0), in 
the limit T p <C u c , pu{t) is the solution of (fSTj) . Such 
a ^-dependence can be measured experimentally and is 
usually an exponential function with initial value and 
final value given by the equilibrium distribution. Accord- 
ing to (|5"Tj) , the initial slope of pn(t) gives the transition 
rate: T_ = pn(0). Likewise, if the system is initial- 
ized in state |1), one can extract T + in a similar way. 
Plotting the resulting transition rates versus bias e, one 
obtains the tunneling resonant peaks. By fitting the ex- 
perimental data to the shifted Gaussian line-shapes ([55]) 
the parameters e p o, W, and T p can be extracted and from 
(f?7)) . A can be obtained. Such a procedure, performed in 
Ref. |3o, successfully confirmed our theory especially the 
relation between W and e p o- 

If the transition rate T p becomes comparable to the 
environment's characteristic energy u) c , the i-local equa- 
tion ([ST]) will not be adequate to describe the evolution 
of the system. However, if the nonlocality effect is small, 
one can still use the same equation but with T± defined 
by ([53)l. hence In such a case, the peak will not 

be symmetric around its center, with an asymmetry that 
dependens on A. Experimental observation of such an 



asymmetry is an indication of time-delayed response of 
the environment and may provide information about to c . 
It should be reminded that a presence of high frequency 
modes in S(uj) may also lead to deviation from a symmet- 
ric Gaussian MKT peak but such an effect is independent 
of A hence could be easily distinguished from the above 
nonlocal effects. 

Another interesting type of experiment is the Landau- 
Zener transition in which e is a linear function of time 
during the evolution. For that type of evolution, again in 
the T p <C ll> c regime, one can still use (|51[) but with a time 
dependent e. Such a procedure was proved successful in 
providing accurate description of the experimental data 
for flux qubits in Ref. [40l . 

It should be mentioned that the tunneling rate A in our 
formalism may not be independent of e as assumed here. 
In practice, as the double-well potential is tilted, it not 
only affects the relative positions of the energy levels in 
the two wells but also affects the matrix elements between 
them. Such dependence is weak for a small bias, but 
as e becomes large the effect of modulation of A might 
become visible. 



B. Tunneling to or between excited states 

If the energy tilt is large enough so that the ground 
state of the initial well becomes in resonance with an ex- 
cited state of the opposite well, tunneling to the excited 
state can occur. Alternatively, one may initialize the 
system in an excited state in the initial well, via e.g., mi- 
crowave excitation, and make the system tunnel between 
two excited states. It is therefore important to under- 
stand how such a tunneling can be described within the 
present theory. One can generalize the arguments of the 
previous section to calculate the tunneling rate. In this 
case, we need to add intrawell relaxations to the picture. 

In Ref. H3, it was shown that the tunneling rate from 
state \i) in the left well to state \j) in the right well is 
given by 



A 2 - 

r«(e) = -f 



dt e i(e-ep)t-7«|t|-3 w2 t 2 



(58) 



where e is the bias energy with respect to the resonance 
point between \i) and \j), Ay is the tunneling amplitude 
between the two states, and 7y = (7i+7j)/2) with 7$ be- 
ing the intrawell relaxation rates corresponding to state 
If one of the states is the ground state in its own well, 
then its corresponding intrawell relaxation rate is zero. 
The transition rate becomes a convolution of Lorentzian 
and Gaussian functions: 



r«(c) 



A?-7ij 



8irW J- 
A 2 - 

Re 

W 



de>- 



2 /2W 2 



(e - e') 2 + 7, 



y?. 



e±e 



P + ijij 



V2W 



(59) 
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where 



w(x) — e x '[l — erf(— ix)] = 



2e 



e~* alt (60) 



is the complex error function. In the limit 7^ — > 0, the 
shifted Gaussian line-shape is recovered. In the op- 
posite limit, jij 3> W , the peak becomes a Lorentzian 
with width 7.^ . 



C. Multi-channel tunneling 

So far, we have investigated the dynamics of a defi- 
nite single tunneling event between the wells. However, 
as the system's temperature increases, one should expect 
the increase of probability of thermal occupation of the 
excited states of each well. Under such conditions, it 
becomes unknown what single tunneling event will take 
place. Consequently, when predicting the effective tun- 
neling rate between wells, one has to take into account 
the statistics of occupation of excited states and their re- 
spective tunneling probabilities to the opposite well, in 
an ensemble average. The net of this thermally assisted 
dynamics is a multi-channel tunneling, which leads to an 
increase of the measured tunneling rate. As we shall see, 
due to the fast increase of the tunneling amplitude Ay 
between excited states \i) and|j), T does not need to be 
too large for this process to become non-negligible. For 
simplicity we consider zero bias (e = 0) situation in which 
the two potential wells are in resonance. 

Let A„ and r± denote the tunneling amplitude and 
transition rate between the n-th energy levels in the op- 
posite wells, and r± the total transition rates between the 
wells. In thermal equilibrium, the occupation probabil- 
ity of the n-th state is given by Boltzmann distribution: 
P n = e- E -l T /Y ii e-- Bi,T . Therefore 



r ± (e) = £p„ii( e ) 



(61) 



At small enough T, one can assume P n w e~ En °/ T (for 
n > 0), where E n o = E n — Eq is the relative energy of 
state |n) compared to the ground state (n = 0). If %j -C 
W for the low-lying energy levels, we may neglect 7^ in 
(|55|) and all r„ will have the same Gaussian functional 
form, leading to 



r-(e) = £< 



8 W 



W 



(62) 



where 



A 2 

n>l 



E n0 /T 



1/2 



(63) 



Therefore, the net contribution from tunneling events in- 
volving excited states can be seen as a renormalization 
of the tunneling amplitude between wells. Since usually 
A„ Ao, such contribution becomes important even at 
temperatures much smaller than the plasma frequency 
uj. p = Eiq. The crossover temperature T co can be ob- 
tained by requiring (Ai/Ao) 2 e~ Wl, / T ~ 1, such that the 
contribution from the first excited state becomes impor- 
tant: 



T r , 



21n(Ai/A ) 



(64) 



Typically Ai is a few orders of magnitude larger than A 
and therefore T co can be an order of magnitude smaller 
than Up. High frequency modes of environment may also 
renormalize the tunneling amplitude^ resulting in a In- 
dependent A e f f even at T < T co . Such a T-dependence is 
typically much weaker and a crossover to the exponential 
dependence in (|63|) should be observable. 



VII. CONCLUSIONS 

We have shown a systematic procedure to determine 
the evolution of a two-state system in the regime of inco- 
herent quantum dynamics. Considering a second order 
perturbation theory in the system bare tunneling rate A, 
and a Gaussian distribution for the environment fluctu- 
ations, we have determined the short time evolution of 
the system reduced density matrix elements. 

Under the assumption of high integrated noise W, i.e., 
a system-bath strong coupling regime, we verify that, 
indeed, dephasing process takes place early in the system 
evolution, which sets 1/W as the smallest time scale of 
the evolution, justifying the claim of having a system 
with incoherent dynamics. 

As for the system populations, we have seen that, in 
general, one should expect complex non-Markovian dy- 
namics. We were able to clearly demonstrate how the 
non-Markovian evolution can be related to the time re- 
sponse of the environment, r^. Indeed, we have veri- 
fied that for time scales t ^> tr, the system follows the 
detailed balance dynamics. On the other hand, if the 
environment response is very slow, i.e., t <C tr, the sys- 
tem sees a static (classical) noise source. In addition, 
by investigating the equation of motion for the reduced 
density matrix, we have demonstrated how one can sim- 
plify the non-Markovian effects by introducing modified 
transition rates for the dynamical equations. 

Finally, we have inspected the intra- and interwell 
transition possibilities inside a double- well potential, and 
quantified how the multi-channel process can lead to an 
enhancement of the system tunneling. We have deter- 
mined the condition for this process to take place, and 
estimated the crossover temperature which can be an or- 
der of magnitude smaller than the system plasma fre- 
quency U! p . 
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Some of the predictions of our theory have already been 
confirmed experimentally 38 ' 40 More experiments, how- 
ever, are necessary especially to confirm our description 
of non-Markovian dynamics. A simple measure of the 
asymmetry of the MKT peak in large A regime could be 
indicative of nonlocality in t. As described in Sec. V, such 
an asymmetry should be A dependent and should dis- 
appear at small A. A A-independent asymmetry could 
result from high frequency components of the environ- 
mental noise that make small lot expansion in (|21[) fail. 
Moreover, a T-dependent measure of the tunneling rates 
can reveal the renormalization of the effective tunnel- 



ing amplitude A due to high frequency noise and the 
crossover temperature T co to the multichannel tunneling 
regime as described in section VI. 
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